library(foreign)
library(survival)
library(lme4)
library(arm)
library(texreg)

cl   <- function(dat,fm, cluster){
           require(sandwich, quietly = TRUE)
           require(lmtest, quietly = TRUE)
           M <- length(unique(cluster))
           N <- length(cluster)
           K <- fm$rank
           dfc <- (M/(M-1))*((N-1)/(N-K))
           uj  <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
           vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)
           coeftest(fm, vcovCL) }
           

kenya <- read.csv("Replication_Data_Kenya.csv")

moi <- kenya[kenya$moi == 1,]
kibaki <- kenya[kenya$kibaki == 1,]

###############################################
########KENYA: TABLE 1 DESCIPTIVE STATS########
###############################################
#Table 1 compiled using the mean() and sd() functions for each regression variable from the kenya dataset

###############################################
##################KENYA: POSTINGS##############
###############################################
native.moi <- glm(non_coethnic_dc ~  non_pres_coethnic_district + joined_normalized +  lpop + larea  , family = "binomial", dat = moi) 
summary(native.moi)
cl.native.moi <- cl(moi, native.moi, moi$district_id)
cl.native.moi

native.moi.lagged.dv  <- glm(non_coethnic_dc ~  non_pres_coethnic_district + joined_normalized +  lpop + larea + dc_lagged , family = "binomial", dat = moi) 
summary(native.moi.lagged.dv)
cl.native.moi.lagged.dv <- cl(moi, native.moi.lagged.dv, moi$district_id)
cl.native.moi.lagged.dv

native.kibaki <- glm(non_coethnic_dc ~  non_pres_coethnic_district + kibaki_appointee +  lpop  + larea  , family = "binomial", dat = kibaki) 
summary(native.kibaki)
cl.native.kibaki <- cl(kibaki, native.kibaki, kibaki$district)
cl.native.kibaki

native.kibaki.lagged.dv  <- glm(non_coethnic_dc ~  non_pres_coethnic_district + kibaki_appointee +  lpop   + larea + dc_lagged , family = "binomial", dat = kibaki) 
summary(native.kibaki.lagged.dv)
cl.native.kibaki.lagged.dv <- cl(kibaki, native.kibaki.lagged.dv, kibaki$district)
cl.native.kibaki.lagged.dv

###############################################
##################KENYA: TENURE################
###############################################
tenure.moi <- coxph(Surv(startpost, time_id, shuffle_nocens) ~ non_pres_coethnic_district +  joined_normalized + non_coethnic_dc +  lpop +  larea, data=moi, method="efron") 
summary(tenure.moi)

tenure.moi.frailty <- coxph(Surv(startpost, time_id, shuffle_nocens) ~  non_pres_coethnic_district +  joined_normalized + non_coethnic_dc +  lpop +  larea + frailty(district_id, distribution="gaussian"), data=moi, method="efron")
 summary(tenure.moi.frailty)

tenure.kibaki <- coxph(Surv(startpost, time_id, shuffle_nocens) ~ non_pres_coethnic_district +  kibaki_appointee + non_coethnic_dc +  lpop +  larea, data=kibaki, method="efron") 
summary(tenure.kibaki)

tenure.kibaki.frailty <- coxph(Surv(startpost, time_id, shuffle_nocens) ~  non_pres_coethnic_district +  kibaki_appointee + non_coethnic_dc +  lpop +  larea + frailty(district_id, distribution="gaussian"), data=kibaki, method="efron")
 summary(tenure.kibaki.frailty)

###############################################
#####KENYA: PREDICTED PROBABILITIES############
###############################################
fm <- native.moi	
dat <- moi
cluster <- factor(moi$district_id)

           require(sandwich, quietly = TRUE)
           require(lmtest, quietly = TRUE)
           M <- length(unique(cluster))
           N <- length(cluster)
           K <- fm$rank
           dfc <- (M/(M-1))*((N-1)/(N-K))
           uj  <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
           vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)
           coeftest(fm, vcovCL) 

native.moi.outcomes1 <- matrix(data=NA, ncol=1, nrow=3) 
native.moi.outcomes2 <- matrix(data=NA, ncol=1, nrow=3) 
set.seed(48109)
N <- 1000
betas <- mvrnorm(n=1000,mu=summary(fm)$coefficients[,1], vcovCL)
dim(betas)

xt.native.moi  <- cbind(1,
		1, 
		mean(dat$joined_normalized),
		mean(dat$lpop),
		mean(dat$larea)
		)

xc.native.moi  <- cbind(1,
		0, 
		mean(dat$joined_normalized),
		mean(dat$lpop),
		mean(dat$larea)
		)

tprob<-apply((1/(1+exp(-(as.matrix(xt.native.moi))%*% t(betas)))),1,as.vector)
tprobA <- apply(tprob,1,mean)

cprob<-apply((1/(1+exp(-(as.matrix(xc.native.moi))%*% t(betas)))),1,as.vector)
cprobA <- apply(cprob,1,mean)

diffprob<-tprobA-cprobA

native.moi.outcomes1[1,1]<- mean(tprobA)
native.moi.outcomes1[2,1]<- quantile(tprobA,.025)
native.moi.outcomes1[3,1]<- quantile(tprobA,.975)

native.moi.outcomes2[1,1]<- mean(cprob)
native.moi.outcomes2[2,1]<- quantile(cprob,.025)
native.moi.outcomes2[3,1]<- quantile(cprob,.975)

fm <- native.moi.lagged.dv	
dat <- moi
cluster <- factor(moi$district_id)

           require(sandwich, quietly = TRUE)
           require(lmtest, quietly = TRUE)
           M <- length(unique(cluster))
           N <- length(cluster)
           K <- fm$rank
           dfc <- (M/(M-1))*((N-1)/(N-K))
           uj  <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
           vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)
           coeftest(fm, vcovCL) 

native.moi.lagged.dv.outcomes1 <- matrix(data=NA, ncol=1, nrow=3) 
native.moi.lagged.dv.outcomes2 <- matrix(data=NA, ncol=1, nrow=3) 
set.seed(48109)
N <- 1000
betas <- mvrnorm(n=1000,mu=summary(fm)$coefficients[,1], vcovCL)
dim(betas)

xt.native.moi.lagged.dv  <- cbind(1,
		1, 
		mean(dat$joined_normalized),
		mean(dat$lpop),
		mean(dat$larea),
		mean(dat$dc_lagged)
		)

xc.native.moi.lagged.dv  <- cbind(1,
		0, 
		mean(dat$joined_normalized),
		mean(dat$lpop),
		mean(dat$larea),
		mean(dat$dc_lagged)
		)

tprob<-apply((1/(1+exp(-(as.matrix(xt.native.moi.lagged.dv))%*% t(betas)))),1,as.vector)
tprobA <- apply(tprob,1,mean)

cprob<-apply((1/(1+exp(-(as.matrix(xc.native.moi.lagged.dv))%*% t(betas)))),1,as.vector)
cprobA <- apply(cprob,1,mean)

diffprob<-tprobA-cprobA

native.moi.lagged.dv.outcomes1[1,1]<- mean(tprobA)
native.moi.lagged.dv.outcomes1[2,1]<- quantile(tprobA,.025)
native.moi.lagged.dv.outcomes1[3,1]<- quantile(tprobA,.975)

native.moi.lagged.dv.outcomes2[1,1]<- mean(cprob)
native.moi.lagged.dv.outcomes2[2,1]<- quantile(cprob,.025)
native.moi.lagged.dv.outcomes2[3,1]<- quantile(cprob,.975)



fm <- native.kibaki	
dat <- kibaki
cluster <- factor(kibaki$district_id)

           require(sandwich, quietly = TRUE)
           require(lmtest, quietly = TRUE)
           M <- length(unique(cluster))
           N <- length(cluster)
           K <- fm$rank
           dfc <- (M/(M-1))*((N-1)/(N-K))
           uj  <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
           vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)
           coeftest(fm, vcovCL) 

native.kibaki.outcomes1 <- matrix(data=NA, ncol=1, nrow=3) 
native.kibaki.outcomes2 <- matrix(data=NA, ncol=1, nrow=3) 
set.seed(48109)
N <- 1000
betas <- mvrnorm(n=1000,mu=summary(fm)$coefficients[,1], vcovCL)
dim(betas)

xt.native.kibaki  <- cbind(1,
		1, 
		mean(dat$kibaki_appointee),
		mean(dat$lpop),
		mean(dat$larea)
		)

xc.native.kibaki  <- cbind(1,
		0, 
		mean(dat$kibaki_appointee),
		mean(dat$lpop),
		mean(dat$larea)
		)

tprob<-apply((1/(1+exp(-(as.matrix(xt.native.kibaki))%*% t(betas)))),1,as.vector)
tprobA <- apply(tprob,1,mean)

cprob<-apply((1/(1+exp(-(as.matrix(xc.native.kibaki))%*% t(betas)))),1,as.vector)
cprobA <- apply(cprob,1,mean)

diffprob<-tprobA-cprobA


native.kibaki.outcomes1[1,1]<- mean(tprobA)
native.kibaki.outcomes1[2,1]<- quantile(tprobA,.025)
native.kibaki.outcomes1[3,1]<- quantile(tprobA,.975)

native.kibaki.outcomes2[1,1]<- mean(cprob)
native.kibaki.outcomes2[2,1]<- quantile(cprob,.025)
native.kibaki.outcomes2[3,1]<- quantile(cprob,.975)


fm <- native.kibaki.lagged.dv	
dat <- kibaki
cluster <- factor(kibaki$district_id)

           require(sandwich, quietly = TRUE)
           require(lmtest, quietly = TRUE)
           M <- length(unique(cluster))
           N <- length(cluster)
           K <- fm$rank
           dfc <- (M/(M-1))*((N-1)/(N-K))
           uj  <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
           vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)
           coeftest(fm, vcovCL) 

native.kibaki.lagged.dv.outcomes1 <- matrix(data=NA, ncol=1, nrow=3) 
native.kibaki.lagged.dv.outcomes2 <- matrix(data=NA, ncol=1, nrow=3) 
set.seed(48109)
N <- 1000
betas <- mvrnorm(n=1000,mu=summary(fm)$coefficients[,1], vcovCL)
dim(betas)

xt.native.kibaki.lagged.dv  <- cbind(1,
		1, 
		mean(dat$kibaki_appointee),
		mean(dat$lpop),
		mean(dat$larea),
		mean(dat$dc_lagged)
		)

xc.native.kibaki.lagged.dv  <- cbind(1,
		0, 
		mean(dat$kibaki_appointee),
		mean(dat$lpop),
		mean(dat$larea),
		mean(dat$dc_lagged)
		)

tprob<-apply((1/(1+exp(-(as.matrix(xt.native.kibaki.lagged.dv))%*% t(betas)))),1,as.vector)
tprobA <- apply(tprob,1,mean)

cprob<-apply((1/(1+exp(-(as.matrix(xc.native.kibaki.lagged.dv))%*% t(betas)))),1,as.vector)
cprobA <- apply(cprob,1,mean)

diffprob<-tprobA-cprobA



native.kibaki.lagged.dv.outcomes1[1,1]<- mean(tprobA)
native.kibaki.lagged.dv.outcomes1[2,1]<- quantile(tprobA,.025)
native.kibaki.lagged.dv.outcomes1[3,1]<- quantile(tprobA,.975)

native.kibaki.lagged.dv.outcomes2[1,1]<- mean(cprob)
native.kibaki.lagged.dv.outcomes2[2,1]<- quantile(cprob,.025)
native.kibaki.lagged.dv.outcomes2[3,1]<- quantile(cprob,.975)

###############################################
#####KENYA: SURVIVAL PLOTS############
###############################################
newdat <- data.frame(non_pres_coethnic_district =c(0,1), joined_normalized =rep(mean(moi$joined_normalized),2), non_coethnic_dc =c(0,0),  lpop =rep(mean(moi$lpop),2), larea =rep(mean(moi$lpop),2))
pdf(file="MoiTenureHazard.pdf", height = 6, width = 6)
par(mar=c(4,4,0,1), bg=NA)
plot(survfit(tenure.moi, newdata=newdat), lwd=5, lty=c(1,3), xlab="", ylab="", main="",  axes=F, col=c("Black"))
axis(1, at=seq(from=0, to=34, by=2), labels=c("0","1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16", "17"), las=1, tick=T, cex.axis=1, col="Black", col.ticks="Black", col.axis="Black")
mtext("Years Since Appointment to District", 1, line=3, col="Black", cex=1.3)
axis(2, cex.axis=1.2, col="Black", col.ticks="Black", col.axis="Black", las=2)
mtext("Proportion of DCs Not Terminated", 2, line=3, col="Black", cex=1.3)
text(3.4, 0.95, "Moi's Core Districts", cex=1.2)
text(3, 0.2, "Districts Outside 
Moi's Core", cex=1.2)
dev.off()


tenure.kibaki <- coxph(Surv(startpost_normalizedkibaki, time_id_normalizedkibaki, shuffle_nocens) ~ non_pres_coethnic_district +  kibaki_appointee + non_coethnic_dc +  lpop +  larea, data=kibaki, method="efron") 
summary(tenure.kibaki)

newdat <- data.frame(non_pres_coethnic_district =c(0,1), kibaki_appointee = c(0,0), non_coethnic_dc =c(0,0), lpop =rep(mean(kibaki$lpop),2),  larea =rep(mean(kibaki$larea),2))
pdf(file="KibakiTenureHazard.pdf", height = 6, width = 6)
par(mar=c(4,4,0,1), bg=NA)
plot(survfit(tenure.kibaki, newdata=newdat), lwd=5, lty=c(1,3), xlab="", ylab="", main="",  axes=F, col=c("Black"))
axis(1, at=seq(from=0, to=34, by=2), labels=c("0","1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16", "17"), las=1, tick=T, cex.axis=1, col="Black", col.ticks="Black", col.axis="Black")
mtext("Years Since Appointment to District", 1, line=3, col="Black", cex=1.3)
axis(2, cex.axis=1.2, col="Black", col.ticks="Black", col.axis="Black", las=2)
mtext("Proportion of DCs Not Terminated", 2, line=3, col="Black", cex=1.3)
text(3.4, 0.95, "Kibaki's Core Districts", cex=1.2)
text(3, 0.2, "Districts Outside 
Kibaki's Core", cex=1.2)
dev.off()

###############################################
##################KENYA: SI - ACLED############
###############################################
native.moi.acled <- glm(non_coethnic_dc ~  non_pres_coethnic_district*acled + joined_normalized +  lpop + larea, family = "binomial", dat = moi) 
summary(native.moi.acled)
cl.native.moi.acled <- cl(moi, native.moi.acled, moi$district_id)
cl.native.moi.acled

native.moi.lagged.dv.acled  <- glm(non_coethnic_dc ~  non_pres_coethnic_district*acled + joined_normalized +  lpop + larea + dc_lagged, family = "binomial", dat = moi) 
summary(native.moi.lagged.dv.acled)
cl.native.moi.lagged.dv.acled <- cl(moi, native.moi.lagged.dv.acled, moi$district_id)
cl.native.moi.lagged.dv.acled

native.kibaki.acled <- glm(non_coethnic_dc ~  non_pres_coethnic_district*acled + kibaki_appointee + lpop + larea, family = "binomial", dat = kibaki) 
summary(native.kibaki)
cl.native.kibaki.acled <- cl(kibaki, native.kibaki.acled, kibaki$district)
cl.native.kibaki.acled

native.kibaki.lagged.dv.acled  <- glm(non_coethnic_dc ~ non_pres_coethnic_district*acled + kibaki_appointee + lpop + larea + dc_lagged, family = "binomial", dat = kibaki) 
summary(native.kibaki.lagged.dv.acled)
cl.native.kibaki.lagged.dv.acled <- cl(kibaki, native.kibaki.lagged.dv.acled, kibaki$district)
cl.native.kibaki.lagged.dv.acled

tenure.moi.acled <- coxph(Surv(startpost, time_id, shuffle_nocens) ~ non_pres_coethnic_district * acled +  joined_normalized + non_coethnic_dc + lpop + larea, data=moi, method="efron") 
summary(tenure.moi.acled)

tenure.moi.frailty.acled <- coxph(Surv(startpost, time_id, shuffle_nocens) ~  non_pres_coethnic_district * acled + joined_normalized + non_coethnic_dc + lpop +  larea + frailty(district_id, distribution="gaussian"), data=moi, method="efron")
summary(tenure.moi.frailty.acled)

tenure.kibaki.acled <- coxph(Surv(startpost, time_id, shuffle_nocens) ~ non_pres_coethnic_district * acled +  kibaki_appointee + non_coethnic_dc + lpop + larea, data=kibaki, method="efron") 
summary(tenure.kibaki.acled)

tenure.kibaki.frailty.acled <- coxph(Surv(startpost, time_id, shuffle_nocens) ~  non_pres_coethnic_district * acled +  kibaki_appointee + non_coethnic_dc +  lpop +  larea + frailty(district_id, distribution="gaussian"), data=kibaki, method="efron")
 summary(tenure.kibaki.frailty.acled)
 
 
###############################################
##############KENYA: SI - LARGER ETH GROUP#####
###############################################
native.moi.larger.group <- glm(non_coethnic_dc ~  non_pres_coethnic_district + joined_normalized +  lpop + larea  , family = "binomial", dat = moi) 
summary(native.moi.larger.group)
cl.native.moi.larger.group <- cl(moi, native.moi.larger.group, moi$district_id)
cl.native.moi.larger.group

native.moi.lagged.dv.larger.group  <- glm(non_coethnic_dc ~  non_pres_coethnic_district + joined_normalized +  lpop + larea + dc_lagged , family = "binomial", dat = moi) 
summary(native.moi.lagged.dv.larger.group)
cl.native.moi.lagged.dv.larger.group <- cl(moi, native.moi.lagged.dv.larger.group, moi$district_id)
cl.native.moi.lagged.dv.larger.group

native.kibaki.larger.group <- glm(non_coethnic_dc ~  non_pres_coethnic_district + kibaki_appointee +  lpop  + larea  , family = "binomial", dat = kibaki) 
summary(native.kibaki.larger.group)
cl.native.kibaki.larger.group <- cl(kibaki, native.kibaki.larger.group, kibaki$district)
cl.native.kibaki.larger.group

native.kibaki.lagged.dv.larger.group  <- glm(non_coethnic_dc ~  non_pres_coethnic_district + kibaki_appointee +  lpop   + larea + dc_lagged , family = "binomial", dat = kibaki) 
summary(native.kibaki.lagged.dv.larger.group)
cl.native.kibaki.lagged.dv.larger.group <- cl(kibaki, native.kibaki.lagged.dv.larger.group, kibaki$district)
cl.native.kibaki.lagged.dv.larger.group

tenure.moi.larger.group <- coxph(Surv(startpost, time_id, shuffle_nocens) ~ non_pres_coethnic_district +  joined_normalized + non_coethnic_dc +  lpop +  larea, data=moi, method="efron") 
summary(tenure.moi.larger.group)

tenure.moi.frailty.larger.group <- coxph(Surv(startpost, time_id, shuffle_nocens) ~  non_pres_coethnic_district +  joined_normalized + non_coethnic_dc +  lpop +  larea + frailty(district_id, distribution="gaussian"), data=moi, method="efron")
 summary(tenure.moi.frailty.larger.group)

tenure.kibaki.larger.group <- coxph(Surv(startpost, time_id, shuffle_nocens) ~ non_pres_coethnic_district +  kibaki_appointee + non_coethnic_dc +  lpop +  larea, data=kibaki, method="efron") 
summary(tenure.kibaki.larger.group)

tenure.kibaki.frailty.larger.group <- coxph(Surv(startpost, time_id, shuffle_nocens) ~  non_pres_coethnic_district +  kibaki_appointee + non_coethnic_dc +  lpop +  larea + frailty(district_id, distribution="gaussian"), data=kibaki, method="efron")
summary(tenure.kibaki.frailty.larger.group)
 


